/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id: AliAnalysisTaskESDMuonFilterO.cxx 57113 2012-06-12 08:52:42Z fca $ */

//
// Add the muon tracks to the generic AOD track branch during the 
// filtering of the ESD. 
//
// Authors: R. Arnaldi 5/5/08 and L. Aphecetche January 2011
//
// Note that we :
//   - completely disable all the branches that are not required by (most) the muon analyses,
//     e.g. cascades, v0s, kinks, jets, etc...
//   - filter severely the tracks (keep only muon tracks) and vertices (keep only primary -including
//     pile-up - vertices) branches 
// 
// (see AddFilteredAOD method)
//

#include "AliAnalysisTaskESDMuonFilterO.h"

#include "AliAODDimuon.h"
#include "AliAODEvent.h"
#include "AliAODHandler.h"
#include "AliAODExtension.h"
#include "AliAODMCParticle.h"
#include "AliAODMuonReplicator0.h"
#include "AliAODVertex.h"
#include "AliAnalysisFilter.h"
#include "AliAnalysisManager.h"
#include "AliCodeTimer.h"
#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "AliESDMuonTrack.h"
#include "AliESDVertex.h"
#include "AliESDtrack.h"
#include "AliLog.h"
#include "AliMCEvent.h"
#include "AliMCEventHandler.h"
#include "AliMultiplicity.h"
#include <TChain.h>
#include <TFile.h>
#include <TParticle.h>
#include "AliESDtrackCuts.h"
#include "AliAODTrack.h"
#include "AliAODPid.h"
#include "TF1.h"
#include "AliPIDResponse.h"
#include "AliAOD3LH.h"

#include "AliAnalysisFilter.h"
#include "AliAODRecoDecay.h"
#include "AliAODRecoDecayHF.h"
#include "AliAODRecoDecayHF2Prong.h"

#include "AliVEvent.h"
#include "AliVVertex.h"
#include "AliVTrack.h"
#include "AliVertexerTracks.h"
#include "AliKFVertex.h"
#include "AliAnalysisFilter.h"
#include "AliAnalysisVertexingHF.h"

using std::cout;
using std::endl;
ClassImp(AliAnalysisTaskESDMuonFilterO)
ClassImp(AliAnalysisNonMuonTrackCuts)

////////////////////////////////////////////////////////////////////////

AliAnalysisNonMuonTrackCuts::AliAnalysisNonMuonTrackCuts()
{
  // default ctor 
}

Bool_t AliAnalysisNonMuonTrackCuts::IsSelected(TObject* obj)
{
  // Returns true if the object is a muon track
  AliAODTrack* track = dynamic_cast<AliAODTrack*>(obj);

  ULong_t  status;

  if(track){
    
    status  = (ULong_t)track->GetStatus();

    if(track->GetTPCNcls() > 80 &&
       track->Chi2perNDF() < 5  &&
       track->IsOn(AliAODTrack::kTPCrefit) &&
       track->IsOn(AliAODTrack::kTPCin)    &&
       !track->IsOn(AliAODTrack::kITSpureSA))
      {
	return kTRUE;
      }
  } 
  
  else 
    return kFALSE;
  

}

AliAnalysisNonPrimaryVertices::AliAnalysisNonPrimaryVertices()
{
  // default ctor   
}

Bool_t AliAnalysisNonPrimaryVertices::IsSelected(TObject* obj)
{
  // Returns true if the object is a primary vertex
  
  AliAODVertex* vertex = dynamic_cast<AliAODVertex*>(obj);
  if (vertex)
    {
      if ( vertex->GetType() == AliAODVertex::kPrimary     ||
	   vertex->GetType() == AliAODVertex::kMainSPD     ||
	   vertex->GetType() == AliAODVertex::kPileupSPD   ||
	   vertex->GetType() == AliAODVertex::kPileupTracks||
	   vertex->GetType() == AliAODVertex::kMainTPC )
	{
	  return kTRUE;
	}
    }
  
  //  enum AODVtx_t {kUndef=-1, kPrimary, kKink, kV0, kCascade, kMulti, kMainSPD, kPileupSPD, kPileupTracks,kMainTPC};

  return kFALSE;
  
}

AliAnalysisTaskESDMuonFilterO::AliAnalysisTaskESDMuonFilterO(Bool_t onlyMuon, Bool_t keepAllEvents, Int_t mcMode, Int_t nsigmaTrk1,Int_t nsigmaTrk2, Int_t partType1,Int_t partType2):
  AliAnalysisTaskSE(),
  fTrackFilter(0x0),
  fEnableMuonAOD(kTRUE),
  fEnableDimuonAOD(kTRUE),
  fOnlyMuon(onlyMuon),
  fKeepAllEvents(keepAllEvents),
  fMCMode(mcMode),
  fnSigmaTrk1(nsigmaTrk1),
  fnSigmaTrk2(nsigmaTrk2),
  fpartType1(partType1),
  fpartType2(partType2),
  fPIDResponse(0),
  fSecVtxWithKF(kFALSE),
  fVertexerTracks(0x0),
  fV1(0x0),
  fAODMapSize(0),
  fAODMap(0),
  fBzkG(0.),
  fCosMin(),
  fDCAtracksMin(),
  fRmax(),
  fRmin(),
  fDNmin(),
  fDPmin()
 
{
  // Default constructor
}

AliAnalysisTaskESDMuonFilterO::AliAnalysisTaskESDMuonFilterO(const char* name, Bool_t onlyMuon, Bool_t keepAllEvents, Int_t mcMode, Int_t nsigmaTrk1,Int_t nsigmaTrk2, Int_t partType1,Int_t partType2):
  AliAnalysisTaskSE(name),
  fTrackFilter(0x0),
  fEnableMuonAOD(kTRUE),
  fEnableDimuonAOD(kTRUE),
  fOnlyMuon(onlyMuon),
  fKeepAllEvents(keepAllEvents),
  fMCMode(mcMode),
  fnSigmaTrk1(nsigmaTrk1),
  fnSigmaTrk2(nsigmaTrk2),
  fpartType1(partType1),
  fpartType2(partType2),
  murep(0),
  fPIDResponse(0),
  fSecVtxWithKF(kFALSE),
  fVertexerTracks(0x0),
  fV1(0x0),
  fAODMapSize(0),
  fAODMap(0),
  fBzkG(0.),
  fCosMin(),
  fDCAtracksMin(),
  fRmax(),
  fRmin(),
  fDNmin(),
  fDPmin()
{
  // Constructor
  
  DefineOutput(1,TList::Class()); // analysis cuts
}

//______________________________________________________________________________
void AliAnalysisTaskESDMuonFilterO::UserCreateOutputObjects()
{
  

  //-----------------------------------------------
  // Particle Identification Setup (new PID object)
  //-----------------------------------------------
  
  AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
  AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
  fPIDResponse = inputHandler->GetPIDResponse();
  
  //--------------------------
  list = new TList();
  list->SetOwner();  // IMPORTANT!

  AODtree = new TTree("AODtree","AODtree");

  fHeader = new TClonesArray("AliAODHeader",30);
  fHeader->SetName("AODHeader");    
 
  fTracks = new TClonesArray("AliAODRecoDecayHF2Prong",30);
  fTracks->SetName("SecondaryVertices");    
  
  fVertices = new TClonesArray("AliAODVertex",2);
  fVertices->SetName("vertices");    
    
  fDimuons = new TClonesArray("AliAODTrack",30);
  fDimuons->SetName("Nuclei");
  
  fVZERO = new TClonesArray("AliAODTrack",30);
  fVZERO->SetName("DaughterTracks");


  AODtree->Branch("AODHeader",  &fHeader);
  AODtree->Branch("SecondaryVertices",  &fTracks);
  AODtree->Branch("vertices",&fVertices);
  AODtree->Branch("Nuclei", &fDimuons);
  AODtree->Branch("DaughterTracks",   &fVZERO);


  // AddAODBranch("TClonesArray", &fTracks  , "cazziemazzi.root");
  // AddAODBranch("TClonesArray", &fVertices, "cazziemazzi.root");
  // AddAODBranch("TClonesArray", &fDimuons , "cazziemazzi.root");
  // AddAODBranch("TClonesArray", &fVZERO   , "cazziemazzi.root");

  list->Add(AODtree);

  //PostData(1,  fListHistCascade);
  PostData(1,  list);
  // Create the output container
  // if (fTrackFilter) OutputTree()->GetUserInfo()->Add(fTrackFilter);
  return;
}

//______________________________________________________________________________
void AliAnalysisTaskESDMuonFilterO::PrintTask(Option_t *option, Int_t indent) const
{
  // Specify how we are configured
  
  AliAnalysisTaskSE::PrintTask(option,indent);
  
  TString spaces(' ',indent+3);
  
  if ( fOnlyMuon ) 
    {
      cout << spaces.Data() << "Keep only muon information " << endl;        
    }
  else 
    {
      cout << spaces.Data() << "Keep all information from standard AOD" << endl;
    }
  
  if ( fKeepAllEvents ) 
    {
      cout << spaces.Data() << "Keep all events, regardless of number of muons" << endl;    
    }
  else 
    {
      cout << spaces.Data() << "Keep only events with at least one muon" << endl;
    }
  
  if ( fMCMode > 0 ) 
    {
      cout << spaces.Data() << "Assuming work on MC data (i.e. will transmit MC branches)" << endl;
    }
}

//______________________________________________________________________________
void AliAnalysisTaskESDMuonFilterO::AddFilteredAOD(const char* aodfilename, const char* title)
{
  // Add an output filtered and replicated aod
  
  AliAODHandler *aodH = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
  if (!aodH) Fatal("UserCreateOutputObjects", "No AOD handler");

  AliAODExtension* ext = aodH->AddFilteredAOD(aodfilename,title);

  if (!ext) return;
  
  if ( fOnlyMuon ) 
    {    
    
      if(!murep)delete murep;

      murep = new AliAODMuonReplicator0("MuonReplicator0",
      					"remove non interesting tracks",
      					new AliAnalysisNonMuonTrackCuts,
      					new AliAnalysisNonPrimaryVertices,
      					fMCMode);
            
      ext->DropUnspecifiedBranches(); // all branches not part of a FilterBranch call (below) will be dropped


      ext->FilterBranch("tracks",murep);    
      ext->FilterBranch("vertices",murep);  
      ext->FilterBranch("dimuons",murep); //per test
      ext->FilterBranch("AliAODVZERO",murep);
      ext->FilterBranch("AliAODTZERO",murep);
      
      cout<<"add filterd aod"<<endl;
      
      if ( fMCMode > 0 ) 
	{
	  // MC branches will be copied (if present), as they are, but only
	  // for events with at least one muon. 
	  // For events w/o muon, mcparticles array will be empty and mcheader will be dummy
	  // (e.g. strlen(GetGeneratorName())==0)
	  
	  ext->FilterBranch("mcparticles",murep);
	  ext->FilterBranch("mcHeader",murep);
	}
    }  
}

//______________________________________________________________________________
void AliAnalysisTaskESDMuonFilterO::Init()
{
  // Initialization
  
  cout<<"Dentro INIT "<<endl;

  // if(fEnableMuonAOD) AddFilteredAOD("AliAOD.Muons.root", "MuonEvents");
  //  if(fEnableDimuonAOD) AddFilteredAOD("AliAOD.Dimuons.root", "DimuonEvents");    
}


//______________________________________________________________________________
void AliAnalysisTaskESDMuonFilterO::UserExec(Option_t */*option*/)
{
  // Execute analysis for current event		       	    
  cout<<"--------------->Enter the user exec<------------------------------"<<endl;

  AliVEvent *event = InputEvent();
  cout<<"Event"<<event<<endl;
  
  // AliAODEvent *lAODevent=(AliAODEvent *)InputEvent();
  AliAODEvent *lAODevent=(AliAODEvent*)InputEvent();

  if (!lAODevent) {
    Printf("ERROR: aod not available");
    return;
  }


  //  if(fEnableMuonAOD) AddFilteredAOD("AliAOD.Muons.root", "MuonEvents");

  

  Long64_t ientry = Entry();
  if(fDebug)printf("Muon Filter: Analysing event # %5d\n", (Int_t) ientry);
  
  
  //Check the PIDresponse
  // if(!fPIDResponse) {
  //   AliError("Cannot get pid response");
  //   return;
  // }

  //***************************************************

  ConvertESDtoAOD(lAODevent);
  
  //*************************************************
  
  if(!murep)delete murep;
  
}

//_____________________________________________________________________

void AliAnalysisTaskESDMuonFilterO::ConvertESDtoAOD(AliAODEvent *source) 

{
  Int_t nheader(0);

  AliAODHandler* handler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
  
  //  AliAODHeader *fHeader0 = source->GetHeader();
  
  new((*fHeader)[nheader++]) AliAODHeader(* source->GetHeader());

  //cout<<"Header :"<<fHeader<<endl;
  // if ( handler /*&& fEnableMuonAOD && (contHe>0)*/) 
  //   {
     
  //     cout<<"\n\n\n------------->QUIIIIIIIIIIIIII: "<<endl;
  //     AliAODExtension *extMuons = handler->GetFilteredAOD("AliAOD.Muons.root");
      
  //     if ( extMuons ){ 
  // 	extMuons->SelectEvent();
  // 	extMuons->Print();
  //     }
  //   }
  
  

  Int_t ntracks(0);
  Int_t ndimuons(0);
  Int_t ntracksD(0);

  Int_t input(0);
  Double_t xdummy,ydummy;

  AliAODRecoDecayHF2Prong *io2Prong  = 0;

  TObjArray *twoTrackArray    = new TObjArray(2);
  Double_t dispersion;

  AliAODVertex *vtx = source->GetPrimaryVertex();
  
  // A Set of very loose cut for a weak strange decay
  
  fCosMin       = 0.97;
  fDCAtracksMin = 1;
  fRmax         = 200.;
  fRmin         = 0.1;
  fDNmin        = 0.05;
  fDPmin        = 0.05;

  //----------------------------------------------------------

  //  Int_t nindices=0;
  UShort_t *indices = 0;
  const Int_t entries = source->GetNumberOfTracks();
  //  AliCentrality* cent;

  Double_t pos[3],cov[6];
  vtx->GetXYZ(pos);
  vtx->GetCovarianceMatrix(cov);
  fV1 = new AliESDVertex(pos,cov,100.,100,vtx->GetName());
  cout<<"fV1 pricipal loop: "<<fV1<<endl;
  
  if(entries<=0) return;
  indices = new UShort_t[entries];
  memset(indices,0,sizeof(UShort_t)*entries);
  fAODMapSize = 100000;
  fAODMap = new Int_t[fAODMapSize];
  memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
  //  cent=((AliAODEvent&)source)->GetCentrality();
  
  //-------------------------------------------------------------

  cout<<"1"<<endl;

  AliAODRecoDecayHF   *rd = 0;

  if(vtx->GetNContributors()<1) {
    
    // SPD vertex cut
    vtx =source->GetPrimaryVertexSPD();
    
    if(vtx->GetNContributors()<1) {
      Info("AliAnalysisTaskHelium3Pi","No good vertex, skip event");
      return; // NO GOOD VERTEX, SKIP EVENT 
    }
  }
  
  Double_t xPrimaryVertex=0.,yPrimaryVertex=0.,zPrimaryVertex=0.;
  xPrimaryVertex=vtx->GetX();
  yPrimaryVertex=vtx->GetY();
  zPrimaryVertex=vtx->GetZ();  
  
  //  Double_t PrimaryVertex[3]={xPrimaryVertex,yPrimaryVertex,zPrimaryVertex};

  //Double_t lMagneticField=source->GetMagneticField();
  fBzkG=source->GetMagneticField();
  fVertexerTracks=new AliVertexerTracks(fBzkG);


  
  TF1 *foPion=new TF1("foPion", "[0]*([1]*TMath::Power(TMath::Sqrt(1 + (x/0.13957)*(x/0.13957))/(x/0.13957) , [3]) - 1 - TMath::Power(TMath::Sqrt(1 + (x/0.13957)*(x/0.13957))/(x/0.13957) , [3])*TMath::Log([2] + 1/TMath::Power((x/0.13957), [4])))",0.01,20);
  foPion->SetParameters(4.1,8.98482806165147636e+00,1.54000000000000005e-05,2.30445734159456084e+00,2.25624744086878559e+00);
  
  Double_t TrackNumber = source->GetNumberOfTracks();
  Int_t    countPionBuoni=0;
  
  cout<<"TrackNumber: "<<TrackNumber<<endl;

  //Tracks arrays
  
  TArrayI Track0(TrackNumber);        //Pions                                                                          
  Int_t nTrack0=0;
  
  TArrayI Track1(TrackNumber);        //Helium3
  Int_t nTrack1=0;
  
  

  for(Int_t j=0; j<TrackNumber; j++){

    //    cout<<"Inside loop tracks"<<endl;
    
    AliVTrack *track = (AliVTrack*)source->GetTrack(j);
    
    AliAODTrack *aodtrack =(AliAODTrack*)track;

    //-----------------------------------------------------------
    //Track cuts 
    if(aodtrack->GetTPCNcls() < 80 )continue;
    if(aodtrack->Chi2perNDF() > 5 )continue;

    // cout<<"Filter Map: "<<aodtrack->GetFilterMap()<<endl;

    //    if (!aodtrack->TestFilterMask(BIT(4))) continue;
    // if(aodtrack->TestFilterMask(BIT(4))!=0)
    //   cout<<"!!!!!!!!!!!!!!!!!!!!!!!!!!! Filter bit mask: !!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
    
    if (!aodtrack->IsOn(AliAODTrack::kTPCrefit)) continue;
    if (!aodtrack->IsOn(AliAODTrack::kTPCin)) continue;
    if (aodtrack->IsOn(AliAODTrack::kITSpureSA)) continue;

    //---------------------------------------------------------------
     
    Double_t mom = aodtrack->GetDetPid()->GetTPCmomentum();
        
    if(mom<0.150)continue;
    
    // Double_t nSigmaNegPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodtrack,(AliPID::EParticleType)fpartType1)); //to be changed
    
    Double_t nSigmaNegPion = TMath::Abs((aodtrack->GetTPCsignal() - foPion->Eval(mom))/foPion->Eval(mom))/0.07;
    
    if ( (nSigmaNegPion <= fnSigmaTrk1) /*&& aodtrack->Charge()==-1*/){
      Track0[nTrack0++]=j;
    }
    
    // Double_t nSigmaNuclei = fPIDResponse->NumberOfSigmasTPC(aodtrack,(AliPID::EParticleType)fpartType2);  //To be changed
    // if(nSigmaNuclei>-fnSigmaTrk2 && aodtrack->GetTPCsignal()<1000 && mom>0.2){ 
    
    Double_t triggerDeDx = 4*AliExternalTrackParam::BetheBlochAleph((mom*2)/(0.938*3),1.0288,31.9806,5.04114e-11,2.13096,2.38541);
    
    if(aodtrack->GetTPCsignal() > triggerDeDx && aodtrack->GetTPCsignal()<5000 && mom>0.2 /*&& aodtrack->Charge()==1*/){
      Track1[nTrack1++]=j;
      new((*fDimuons)[ndimuons++]) AliAODTrack(*aodtrack);
   
      // cout<<"-------------------------------------------------------------------->fDimuons"<<fDimuons<<endl;
      // cout<<"AOD track: "<<aodtrack<<endl;
      //      fDimuons->Delete();
    }
  
  } 


  
  //cout<<"Delete AodTrack..."<<endl;
  //  delete aodtrack;
  //  aodtrack->Delete();

  Track0.Set(nTrack0);
  Track1.Set(nTrack1);
  
  cout<<"Track loops..."<<endl;
  
  AliAODTrack *postrackAOD = 0;
  AliAODTrack *negtrackAOD = 0;
 
  AliESDtrack *postrack = 0;
  AliESDtrack *negtrack = 0;

  cout <<"nTrack1 "<< nTrack1<<endl;

  for (Int_t i=0; i<nTrack1; i++){                            //! He Tracks Loop

    cout<<"Inside track loop"<<endl;
    
    Int_t Track1idx=Track1[i];

    AliVTrack *trackpos = (AliVTrack*)source->GetTrack(Track1idx);
    postrackAOD =(AliAODTrack*)trackpos;
    postrack = new AliESDtrack(trackpos);
    
    // postrackAOD=source->GetTrack(Track1idx); //Questo

    for (Int_t k=0; k <nTrack0 ; k++) {                           //! Pion Tracks Loop

      Int_t Track0idx=Track0[k];

      AliVTrack *trackneg = (AliVTrack*)source->GetTrack(Track0idx);
      negtrackAOD =(AliAODTrack*)trackneg;
      negtrack = new AliESDtrack(trackneg); //<---questo lo turba proprio parecchio...
      
      //      negtrackAOD=source->GetTrack(Track0idx); //Questo
      
      // negtrack = new AliESDtrack((AliVTrack*)negtrackAOD); //Questo
      // postrack = new AliESDtrack((AliVTrack*)postrackAOD); //Questo
   
      twoTrackArray->AddAt(negtrack,0);
      twoTrackArray->AddAt(postrack,1);
      
      Double_t dcap1n1 = postrack->GetDCA(negtrack,fBzkG,xdummy,ydummy);

      Double_t dcap1toPV = TMath::Abs(postrack->GetD(xPrimaryVertex, yPrimaryVertex,fBzkG));
      Double_t dcan1toPV = TMath::Abs(negtrack->GetD(xPrimaryVertex, yPrimaryVertex,fBzkG));

      if(dcap1n1>fDCAtracksMin)continue;
      if((xdummy+ydummy)>fRmax )continue;
      if((xdummy+ydummy)< fRmin)continue;
      
      if ( dcan1toPV < fDNmin)               
        if ( dcap1toPV < fDPmin) continue;   
      
      cout<<"dcap1n1: "<<dcap1n1<<endl;
 
      Int_t charge[2];
      charge[0]=negtrack->Charge();
      charge[1]=postrack->Charge();
    
      cout<<"Charge neg: "<<negtrack->Charge()<<endl;
      cout<<"Charge pos: "<<postrack->Charge()<<endl;
      
      AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray,dispersion,kTRUE);
      
      if(!vertexp1n1) {

	twoTrackArray->Clear();
	delete negtrack;
	negtrack=NULL; 

      	continue; 
      }

      cout<<"vertexp1n1:  "<<vertexp1n1<<endl;

      io2Prong = Make2Prong(twoTrackArray,source,vertexp1n1,dcap1n1);
  
      if(io2Prong->CosPointingAngle()<fCosMin)continue;

      cout<<"pointing angle "<<io2Prong->CosPointingAngle()<<endl;
      
      rd =  new((*fTracks)[ntracks++]) AliAODRecoDecayHF2Prong(*io2Prong);
      
      
      rd->SetSecondaryVtx(vertexp1n1);
      vertexp1n1->SetParent(rd);
      //      AddRefs(vertexp1n1,rd,source,twoTrackArray);
      
      //vertexp1n1->Delete();  delete vertexp1n1;
	
      new((*fVZERO)[ntracksD++]) AliAODTrack(*negtrackAOD);
      new((*fVZERO)[ntracksD++]) AliAODTrack(*postrackAOD);
      //  fVZERO->Delete();
     
     
      // delete vertexp1n1;
      // vertexp1n1= NULL;
      // continue;
    }

    // delete postrack; 
    // postrack = NULL;
    cout<<"Header :"<<fHeader<<endl;
    cout<<"fVZERO:   "<<    fVZERO<<endl;
    cout<<"fDimuons: "<<    fDimuons<<endl;
    cout<<fHeader   <<endl;
    cout<<fTracks   <<endl;
    cout<<fVertices <<endl;
    cout<<fDimuons  <<endl;
    cout<<fVZERO    <<endl;

   
    AODtree->Fill();
    fHeader   ->Delete();
    fTracks   ->Delete();
    fVertices ->Delete();
    fDimuons  ->Delete();
    fVZERO    ->Delete();
    fHeader   ->Delete();
    // delete negtrack; 
    // negtrack = NULL;
    // delete postrack; 
    // postrack = NULL;
 

  }
  
  //----------------------------------------------------------
  
  assert(fVertices!=0x0);
  fVertices->Clear("C");
  TIter nextV(source->GetVertices());
  AliAODVertex* v;
  Int_t nvertices(0);
  
  while ( ( v = static_cast<AliAODVertex*>(nextV()) ) )
    {
      // if ( !fVertexCut || fVertexCut->IsSelected(v) ) 
      // 	{
	  AliAODVertex* tmp = v->CloneWithoutRefs();
	  AliAODVertex* copiedVertex = new((*fVertices)[nvertices++]) AliAODVertex(*tmp);
	  
	  // to insure the main vertex retains the ncontributors information
	  // (which is otherwise computed dynamically from
	  // references to tracks, which we do not keep in muon aods...)
	  // we set it here
	  
	  copiedVertex->SetNContributors(v->GetNContributors()); 
	  
	  fVertices->Delete();
	  //	  delete copiedVertex;
	  delete tmp;
	  //	}
    }
  
  printf("....Ho eseguito il muon Replicator...\n");
  
  AliDebug(1,Form("input mu tracks=%d tracks=%d vertices=%d ndimuons=%d",
                  input,fTracks->GetEntries(),fVertices->GetEntries(),fDimuons->GetEntries())); 
  
  cout<<"....Sotto debug..."<<endl;
  
  // Finally, deal with MC information, if needed
  
  // if ( fMCMode > 0 )
  //   {
  //     FilterMC(source);      
  //   }
  
  delete foPion;
  foPion = NULL;

  cout<<"Delete..."<<endl;
  cout<<"Delete 1"<<  endl;
  if(io2Prong) {delete io2Prong; io2Prong=NULL;}
  cout<<"Delete 2"<<  endl;
  twoTrackArray->Delete();  delete twoTrackArray;
  cout<<"Delete 3"<<  endl;
  // vtx->Delete();  delete vtx;
  cout<<"Delete 4"<<  endl;
  if(fV1) { delete fV1; fV1=NULL; }
  cout<<"Delete 5"<<  endl;
  if(fAODMap) { delete [] fAODMap; fAODMap=NULL; }
  delete indices;
  cout<<"Delete 6"<<  endl;
  delete fVertexerTracks;

  cout<<"Fine Replicator"<<endl;


}


//-----------------------------------------------------------------------------

AliAODVertex *AliAnalysisTaskESDMuonFilterO::ReconstructSecondaryVertex(TObjArray *trkArray,
								Double_t &dispersion,Bool_t useTRefArray) const
{
  // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
  //AliCodeTimerAuto("",0);

  AliESDVertex *vertexESD = 0;
  AliAODVertex *vertexAOD = 0;

  if(!fSecVtxWithKF) { // AliVertexerTracks

    fVertexerTracks->SetVtxStart(fV1);
    vertexESD = (AliESDVertex*)fVertexerTracks->VertexForSelectedESDTracks(trkArray);

    if(!vertexESD) return vertexAOD;


    Double_t vertRadius2=vertexESD->GetXv()*vertexESD->GetXv()+vertexESD->GetYv()*vertexESD->GetYv();
    if(vertRadius2>200.){
      // vertex outside beam pipe, reject candidate to avoid propagation through material
      delete vertexESD; vertexESD=NULL;
      return vertexAOD;
    }

  } else { // Kalman Filter vertexer (AliKFParticle)

    AliKFParticle::SetField(fBzkG);

    AliKFVertex vertexKF;

    Int_t nTrks = trkArray->GetEntriesFast();
    for(Int_t i=0; i<nTrks; i++) {
      AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
      AliKFParticle daughterKF(*esdTrack,211);
      vertexKF.AddDaughter(daughterKF);
    }
    vertexESD = new AliESDVertex(vertexKF.Parameters(),
				 vertexKF.CovarianceMatrix(),
				 vertexKF.GetChi2(),
				 vertexKF.GetNContributors());

  }
  // convert to AliAODVertex
  Double_t pos[3],cov[6],chi2perNDF;
  vertexESD->GetXYZ(pos); // position
  vertexESD->GetCovMatrix(cov); //covariance matrix
  chi2perNDF = vertexESD->GetChi2toNDF();
  dispersion = vertexESD->GetDispersion();
  delete vertexESD; 
  vertexESD=NULL;

  Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
  vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);

  cout<<"------------------> Reconstruct vertexAOD: "<<vertexAOD<<endl;

  return vertexAOD;

 
}

//-----------------------------------------------------------------------------

 // AliAODRecoDecayHF2Prong* AliAnalysisTaskESDMuonFilterO::Make2Prong(TObjArray *twoTrackArray,const AliAODEvent &evento,
 // 							   AliAODVertex *secVert,Double_t dca)

AliAODRecoDecayHF2Prong* AliAnalysisTaskESDMuonFilterO::Make2Prong(TObjArray *twoTrackArray,AliAODEvent *evento,
							   AliAODVertex *secVert,Double_t dca)

{ 

  cout<<"Inside make 2 prong"<<endl;

  Int_t charge[2];
  charge[0]=-1;
  charge[1]=2;
      
  // From AliAnalysisVertexingHF.cxx Method:Make2Prongs
  
  cout<< charge[0]<<" "<< charge[1]<<endl;
      
  fBzkG = evento->GetMagneticField();
      
  cout<<"fBzkG: "<<fBzkG<<endl;
  cout<<"twoTrackArray: "<<twoTrackArray<<endl;
  Double_t px[2],py[2],pz[2],d0[2],d0err[2];
  AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
  AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
  
  cout<<"negtrack: "<<negtrack<<" postrack: "<<postrack<<endl;
  cout<<"kVeryBig: "<<kVeryBig<<endl;
  cout<<"secVert: "<<secVert<<endl;

  // // propagate tracks to secondary vertex, to compute inv. mass
  postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
  negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
  
  Double_t momentum[3];
  postrack->GetPxPyPz(momentum);
  px[0] = charge[1]*momentum[0]; py[0] = charge[1]*momentum[1]; pz[0] = charge[1]*momentum[2]; 
  cout<< px[0] <<" "<< " "<< py[0] << " "<< pz[0]<<endl;
  //  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2]; 
  negtrack->GetPxPyPz(momentum);
  px[1] = charge[0]*momentum[0]; py[1] = charge[0]*momentum[1]; pz[1] = charge[0]*momentum[2]; 
  cout<< px[1] <<" "<< " "<< py[1] << " "<< pz[1]<<endl;
  //px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2]; 
  

  // primary vertex to be used by this candidate
  AliAODVertex *primVertexAOD  = evento->GetPrimaryVertex();
  cout<<"primVertexAOD "<<primVertexAOD<<endl;
  if(!primVertexAOD) return 0x0;
      
  Double_t d0z0[2],covd0z0[3];
  //  postrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
  d0[0] = d0z0[0];
  d0err[0] = TMath::Sqrt(covd0z0[0]);
  // negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
  d0[1] = d0z0[0];
  d0err[1] = TMath::Sqrt(covd0z0[0]);
  
  // create the object AliAODRecoDecayHF2Prong
  //  AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dcap1n1);
  AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dca);
  cout<<"the2Prong "<<the2Prong<<endl;
  the2Prong->SetOwnPrimaryVtx(primVertexAOD);
  
  UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
  the2Prong->SetProngIDs(2,id);
  //delete primVertexAOD; primVertexAOD=NULL;
 
  if(postrack->Charge()!=0 && negtrack->Charge()!=0) { // don't apply these cuts if it's a Dstar 
      
    //    AddDaughterRefs(secVert,(AliAODEvent&)evento,twoTrackArray);
    AddDaughterRefs(secVert,evento,twoTrackArray);
      
  }
  
  return the2Prong;  

  delete the2Prong;
}

//----------------------------------------------------------------------------
// void AliAnalysisTaskESDMuonFilterO::AddDaughterRefs(AliAODVertex *v,
// 					    const AliAODEvent &event,
// 					    const TObjArray *trkArray) const
void AliAnalysisTaskESDMuonFilterO::AddDaughterRefs(AliAODVertex *v,
					    const AliAODEvent *event,
					    const TObjArray *trkArray) const

{
  // Add the AOD tracks as daughters of the vertex (TRef)
  //AliCodeTimerAuto("",0);
  cout<<"Inside  AddDaughterRefs "<<endl;

  Int_t nDg = v->GetNDaughters();
  
  cout<<"nDg "<<nDg<<endl;

  TObject *dg = 0;
  if(nDg) dg = v->GetDaughter(0);
  cout<<"1 "<<1<<endl;
  if(dg) return; // daughters already added
  cout<<"2 "<<2<<endl;
  Int_t nTrks = trkArray->GetEntriesFast();
  cout<<"nTrks "<<nTrks<<endl;
  AliExternalTrackParam *track = 0;
  AliAODTrack *aodTrack = 0;
  Int_t id;
  
  for(Int_t i=0; i<nTrks; i++) {
    track = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
    cout<<"track inside loop: "<<track<<endl;
    id = (Int_t)track->GetID();
    printf("---> %d\n",id);
    if(id<0) continue; // this track is a AliAODRecoDecay
    cout<<aodTrack<<endl;
    cout<<"1 "<<fAODMap[id]<<endl;
    aodTrack = (AliAODTrack*)event->GetTrack(fAODMap[id]);
    cout<<"2 "<<fAODMap[id]<<endl;
    v->AddDaughter(aodTrack);
    cout<<"End Loop"<<endl;
  }
  cout<<"END AddDaughterRefs"<<endl;
  return;
}
//----------------------------------------------------------------------------
	
void AliAnalysisTaskESDMuonFilterO::AddRefs(AliAODVertex *v,AliAODRecoDecayHF *rd,
				    const AliAODEvent *event,
				    const TObjArray *trkArray) const
// void AliAnalysisTaskESDMuonFilterO::AddRefs(AliAODVertex *v,AliAODRecoDecayHF *rd,
// 				    const AliAODEvent &event,
// 				    const TObjArray *trkArray) const

{
  // Add the AOD tracks as daughters of the vertex (TRef)
  // Also add the references to the primary vertex and to the cuts
  //AliCodeTimerAuto("",0);
  
  AddDaughterRefs(v,event,trkArray);
  //  rd->SetPrimaryVtxRef((AliAODVertex*)event.GetPrimaryVertex());
  rd->SetPrimaryVtxRef(event->GetPrimaryVertex());
  return;
}	

//-----------------------------------------------------------------------------



//------------------------------------------------
void AliAnalysisTaskESDMuonFilterO::Terminate(Option_t */*option*/)
{
  // Terminate analysis
  //
  //  delete murep;
  
  if (fDebug > 1) printf("AnalysisESDfilter: Terminate() \n");
}
